cytofWorkflow之读入FCS文件（一）

前面我们公布了《cytof数据资源介绍（文末有交流群）》，现在就开始正式手把手教学。首先是cytofWorkflow的解析，这个cytofWorkflow基本上涵盖了以FCS文件格式存放的cytof数据处理的方方面面，非常建议初学者从这个开始掌握cytof数据处理！

关于 cytofWorkflow

在bioconductor官网链接是：

• https://bioconductor.org/packages/release/workflows/html/cytofWorkflow.html
• 教程：https://bioconductor.org/packages/release/workflows/vignettes/cytofWorkflow/inst/doc/cytofWorkflow.html
• 开发者主页：http://imlspenticton.uzh.ch/robinson_lab/cytofWorkflow/

实际上跑一下cytofWorkflowbioconductor官网教程就足够了，我这里把他们的教程拓展一下，以一篇发表在nature medicine杂志的文章数据为例子，演示给大家全部的流程细节。

下载FCS文件

前面我们介绍了流式细胞术数据文件标准，FCS文件类型 ，来龙去脉。大家只需要了解最新版即可，是2020年9月3日Spidlen J等人在Cytometry A杂志上提出了的FCS 3.2，该修订版满足了一些新需求和建议，并结合了十年来整个细胞仪领域的进步。FCS 3.2规范完整版可以下载http://flowcyt.sf.net/fcs/fcs32.pdf

当然了，我们如果处理的是公共数据，通常是比较老的版本，不过，只要是FCS 3.0以上就ok。比如发表在nature medicine杂志的文章《Immune profiling of human tumors identifies CD73 as a combinatorial target in glioblastoma》：

• 链接是：https://www.nature.com/articles/s41591-019-0694-x

它的 Data availability 部分就清晰地列出来了cytof数据和单细胞转录组数据存放的地方：

• CyTOF data (Figs. 1, 2a,b, 3c–g and 4c,d) have been deposited with the FlowRepository (FR-FCM-Z2B3).
• scRNA-seq data (Figs. 2c–e and 3a) have been deposited with the Sequence Read Archive with accession number PRJNA588461.

根据链接 http://flowrepository.org/id/FR-FCM-Z2B3 ，可下载到全部的FCS文件。我这里演示其中5个文件，如下：

  3.6M Nov  6 10:20 WT-UnRx 3-G2_CD45pos.fcs
   13M Nov  6 10:27 WT-UnRx-1_CD45pos.fcs
   25M Nov  6 10:27 WT-UnRx-2_CD45pos.fcs
   14M Nov  6 10:27 WT-UnRx-4_CD45pos.fcs
   15M Nov  6 10:20 WT-UnRx-5_CD45pos.fcs

请自行下载哈！

读取FCS文件

只需要把全部的下载成功的FCS文件统一存放在一个文件夹，比如 paper_fcs文件夹，然后使用read.flowSet函数即可读取，示例代码如下：

rm(list = ls())
require(cytofWorkflow) 
p1='paper_fcs/'
fs1=list.files(p1,'*fcs' )
fs1
samp <- read.flowSet(files = fs1,path = p1)
samp

熟悉read.flowSet函数产生的对象

使用read.flowSet函数全部的FCS文件后，会产生一个对象，这里面变量名是 samp 。如下所示：

> samp
A flowSet with 5 experiments.

  column names:
  As75Di Bi209Di Ce140Di Center Cs133Di Dy161Di Dy162Di Dy163Di Dy164Di Er166Di Er167Di Er168Di Er170Di Eu151Di Eu153Di Event_length Gd155Di Gd156Di Gd158Di Gd160Di Ho165Di Ir191Di Ir193Di Lu175Di Nd142Di Nd143Di Nd144Di Nd145Di Nd146Di Nd148Di Nd150Di Offset Os190Di Pb208Di Pd102Di Pd104Di Pd105Di Pd106Di Pd108Di Pd110Di Pr141Di Pt195Di Residual Sm147Di Sm149Di Sm152Di Sm154Di Sn120Di Tb159Di Tm169Di Xe131Di Y89Di Yb171Di Yb172Di Yb173Di Yb174Di Yb176Di beadDist Time
> 

我演示的是读取了5个文件，如下：

  3.6M Nov  6 10:20 WT-UnRx 3-G2_CD45pos.fcs
   13M Nov  6 10:27 WT-UnRx-1_CD45pos.fcs
   25M Nov  6 10:27 WT-UnRx-2_CD45pos.fcs
   14M Nov  6 10:27 WT-UnRx-4_CD45pos.fcs
   15M Nov  6 10:20 WT-UnRx-5_CD45pos.fcs

他们的抗体数量是一致的，所以这个read.flowSet函数没有报错。这个samp对象呢，有5个元素，每个元素也可以独立使用  exprs 函数获取其表达量信息，如下所示：

> exprs(samp[[1]])[1:6, 1:5]
     As75Di    Bi209Di   Ce140Di  Center Cs133Di
[1,]      0 0.00000000 0.0000000 639.176       0
[2,]      0 0.07410723 0.0000000 521.575       0
[3,]      0 0.00000000 0.1077218 716.634       0
[4,]      0 0.00000000 0.0000000 821.956       0
[5,]      0 0.00000000 0.0000000 684.997       0
[6,]      0 0.69866854 0.0000000 786.601       0 

每个元素，都是通过索引获取：

> lapply(1:5, function(x) dim( exprs(samp[[x]]) ))
[[1]]
[1] 15841    59

[[2]]
[1] 58421    59

[[3]]
[1] 109203     59

[[4]]
[1] 62320    59

[[5]]
[1] 67497    59

cytofWorkflow流程后续分析，都是基于这个读取好的对象哦。可以看到，每个样品的细胞数量差异非常大！！！

cytofWorkflow流程示例数据集

其实R包 HDCytoData 就内置了一些cytof数据集哈， 不同数据集，需要不同的函数来下载，所以对网速要求比较高：

library(HDCytoData)
fs <- Bodenmiller_BCR_XL_flowSet()
# 如果网络不好，也可以自行下载
# 然后：loaded into R as a flowSet using read.flowSet() from the flowCore package

这个下载到的fs变量，就就等价于前面的使用read.flowSet函数全部的FCS文件后产生的对象：

> fs
A flowSet with 16 experiments.

  column names:
  Time Cell_length CD3(110:114)Dd CD45(In115)Dd BC1(La139)Dd BC2(Pr141)Dd pNFkB(Nd142)Dd pp38(Nd144)Dd CD4(Nd145)Dd BC3(Nd146)Dd CD20(Sm147)Dd CD33(Nd148)Dd pStat5(Nd150)Dd CD123(Eu151)Dd pAkt(Sm152)Dd pStat1(Eu153)Dd pSHP2(Sm154)Dd pZap70(Gd156)Dd pStat3(Gd158)Dd BC4(Tb159)Dd CD14(Gd160)Dd pSlp76(Dy164)Dd BC5(Ho165)Dd pBtk(Er166)Dd pPlcg2(Er167)Dd pErk(Er168)Dd BC6(Tm169)Dd pLat(Er170)Dd IgM(Yb171)Dd pS6(Yb172)Dd HLA-DR(Yb174)Dd BC7(Lu175)Dd CD7(Yb176)Dd DNA-1(Ir191)Dd DNA-2(Ir193)Dd group_id patient_id sample_id population_id

还有更多其它内置数据集，也可以使用如下所示的函数：

ehub <- ExperimentHub() # create ExperimentHub instance 
ehub <- query(ehub, "HDCytoData") # find HDCytoData datasets 
md <- as.data.frame(mcols(ehub)) # retrieve metadata table

帮助文档写的很清楚：

 

最好的入门方式

如果你也想开启自己的生物信息学数据处理生涯，但是自学起来困难重重，还等什么呢，赶快行动起来吧！参加我们生信技能树官方举办的学习班：

• 数据挖掘学习班第8期（线上直播3周，马拉松式陪伴，带你入门），原价4800的数据挖掘全套课程， 疫情期间半价即可抢购。
• 生信爆款入门-第10期（线上直播4周，马拉松式陪伴，带你入门），原价9600的生信入门全套课程，疫情期间3.3折即可抢购。

生信技能树的粉丝都知道我们有一个全国巡讲的良心学习班，口碑爆棚，生物信息学入门省心省时省力！先看看大家的反馈吧：